A novel formulation of nonlocal electrostatics 
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The accurate modeling of the dielectric properties of water is crucial for many applications in 
physics, computational chemistry and molecular biology. This becomes possible in the framework 
of nonlocal electrostatics, for which we propose a novel formulation allowing for numerical solutions 
for the nontrivial molecular geometries arising in the applications mentioned before. Our approach 
is based on the introduction of a secondary field, ij), which acts as the potential for the rotation free 
part of the dielectric displacement field D. For many relevant models, the dielectric function of the 
medium can be expressed as the Green's function of a local differential operator. In this case, the 
resulting coupled Poisson (-Boltzmann) equations for i/> and the electrostatic potential reduce to 
a system of coupled PDEs. The approach is illustrated by its application to simple geometries. 

PACS numbers: 41.20.Cv, 77.22.Ch 



The theory of continuum electrostatics plays a ma- 
jor role in the determination of solvation free energies of 
atoms, ions, and biomolecules [ij]. Recent progress in its 
applicability to biological systems has been impressive: 
the electrostatic potentials of large biomolecules such as, 
e.g., microtubuli and ribosomes, can be determined 
Unfortunately, the standard continuum approach ulti- 
mately becomes inaccurate when used to determine elec- 
trostatic properties on atomic scales Q, as it is feature- 
less, i.e., the correlation between solvent arrangements 
and the geometrical structure of biomolecular assemblies 
is not taken into account. On the other hand, continuum 
electrostatics is still much more efficient from a computa- 
tional point of view than microscopic simulations based 
on, e.g., molecular dynamics (MD). Therefore, interest 
has risen recently in extensions of the theory of contin- 
uum electrostatics that allow to account for spatial vari- 
ations of the dielectric behaviour of the solvent, in par- 
ticular near boundaries P, Part of the motivation 
for such approaches stems from the field of protein dock- 
ing, where a realistic and efficient modelling of solvent 
properties is essential 

Within the continuum theory of electrodynamics, spa- 
tial dispersion effects can be taken into account in an 
approach called 'nonlocal electrostatics' [ESSSl- I* 
rests on the assumption of a linear relationship between 
the dielectric displacement field and the electric field me- 
diated by a permittivity kernel depending on two spatial 
arguments, 

B{r)^eoJdr'e{r,r')E{r') (1) 
where £(r, r') is the dielectric permittivity tensor. Equiv- 
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alently to eq. one can express the nonlocal relationship 
in terms of the polarization fields |^. 

While the theory of nonlocal electrostatics remains 
firmly embedded within the well-understood framework 
of Maxwell's theory, it introduces a new characteristic 
length scale absent in local electrostatics: the correlation 
length A of the polarization correlations between the sol- 
vent molecules. This length sets the relevant scale for 
the deviation of the dielectric properties of the solvent 
from its bulk value. Thus, nonlocal electrostatics is a se- 
rious candidate for a more realistic description of solvent 
properties, provided it is also computationally tractable. 
It is here where the difficulties arise, however. The the- 
ory of nonlocal electrostatics, discussed in detail below, 
is technically considerably more demanding than local 
electrostatics, since it is usually formulated as a system 
of coupled integro-differential equations. Consequently, 
it has so far only been applied to idealized situations, 
and even then typically after introducing additional ap- 
proximations in order to obtain analytical results [lOlllll| . 
For complex geometries, the solution of the equations by 
numerical methods becomes a formidable task. 

Thus, a reformulation of the equations of nonlocal elec- 
trostatics is needed in order to make the theory applica- 
ble to real-world problems. Here we present a scheme 
which allows to rewrite the theory in terms of a system 
of partial differential equations for the local fields which 
consequently makes it amenable to standard methods of 
numerical analysis. This derivation relies on two assump- 
tions which are typically valid for e.g. the discussion 
of solvation problems for biomolecules: (i) the linearity 
of the relationship between dielectric response and the 
electric field, and (ii) the representation of the dielectric 
function in terms of Green functions of known differential 
operators. We first derive the set of equations, and then 
illustrate how to solve them on (simple) examples. 

In the following we consider the situation given in Fig- 
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FIG. 1: A cavity with local dielectric constant (Q) immersed 
in a medium with spatially varying permittivity (S). The 
interface between the two spaces is denoted by F. 

ure 1. A domain il (which represents a molecule) is em- 
bedded in the solvent E which takes up the whole space 
except rj. The surface of the embedded domain is de- 
noted by r. Within il, the dielectric properties are as- 
sumed to be local, i.e., e = sq, while in the solvent space 

E, eq.([3) holds. On large scales, i.e. when |r — r'| A, 
the dielectric response in the solvent is again local with 
the bulk dielectric constant (with, e.g., es ~ 80 for 
water). In this situation the standard equations of non- 
local electrostatics read as 

A(j)n = —g , reil (2) 

eoV / dr'e(r,r')V>s(r') = 0, r,r'eE (3) 

where g is the density of fixed charges which are as- 
sumed to lie confined within f2\r, which is usually the 
case for biomolecules in solution. This is no restriction 
on the validity of our approach: the existence of surface 
charges on F only slightly modifies the boundary condi- 
tions. Eas. (|2l3() are the Poisson equations for the geom- 
etry of Fig. 1. The primed symbol V' denotes the differ- 
entiation with respect to r'. Without surface charges on 

F, the boundary conditions to cas. H2l3|l for the normal 
(n) and tangential (t) components of the electric and di- 
electric displacement field on the boundary F are given 

by 

Dn.n = Ds,n , En^t = Es^t (4) 

where, by virtue of eq.lPJ, the boundary condition for D 
is also nonlocal. 

In order to step over from the integro-differential to a 
purely differential formulation we introduce, in addition 
to the potential field (j), another potential field tp within 
both compartments and S. Attempts similar in spirit, 
but differing in the implementation, have been discussed 
before in the literature [T^ . again resulting in systems 



of integro-differential equations. In we define the re- 
lations between the potentials 4> a-nd iJj and the physical 
fields as 

En = -^(pn , Do = -Vi/- (5) 
while in E we have 

Ee = -V0s^ (6) 

while the dielectric displacement field in S can be repre- 
sented in terms of a scalar and a vector field [l^ . 

Ds = -eo / dr'eir, r') V'^s = -Vt/^e + V x ■ (7) 

The scalar field ip thus serves as the potential of the rota- 
tion free part D — — Vi/'s of the dielectric displacement 
field D. In our setting 0e is determined by "02 alone, 
and neither ipQ nor (j>Q are affected by ^s. Note that 
can of course be computed from e.g. eq. as soon as 
(^s is known. Since we here are interested only in the 
electrostatic potential (j) and quantities derived from it, 
we will ignore in the following. 

With this definition, the differential equations and 
boundary conditions, eqs.|(51-|7I) can be brought into the 
form 

Aijn = -g , ren, A^^ = , r e E (8) 

and 

VV'sIn = VV-nU , r e F , V(/.s|t = V0n|t , r e F (9) 

In addition, there are now two equations relating the po- 
tential fields (j) and ip, 

eoencpn = ipn , r e O (10) 

and eq.Q, to be fulfilled in E. 

So far, the introduction of t/j has made it possible to 
rewrite the boundary conditions in a completely local 
way. The nonlocal problem now only is with eq.lQ in 
which the nonlocal integral kernel still remains. The lat- 
ter expression can be simplified further under additional 
assumptions on the explicit form of e(r, r'). We now as- 
sume that it can be written in the form 

e(r,r') = e£(5(r-r') +ea(r,r') (11) 

where e — (es — e^)/A^ > 0, and the Green function Q 
solves 

Cg{v,v') = -S{r-r') (12) 

for a differential operator C with constant coefficients. 
Note that ei refers to the value of the dielectric function 
on smallest scales (i.e., r r') and can be related to the 
frequency spectrum of the dielectric function [H . While 
this construction clearly restricts the theory of nonlo- 
cal electrostatics to a certain class of dielectric functions. 
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this restriction is not problematic as it applies to most 
situations of physical interest. 

Under these assumptions the application of the differ- 
ential operator £ to eq.Q yields 



£0 {eiL - e\ V0S -CS/ijj^ . 



(13) 



As the fields and -0 are determined only up to an arbi- 
trary constant, we can drop the gradients on both sides 
of eq. H13|) and are then left with the expression 



(14) 



Eas.l|51- llUll4|) constitute the novel formulation of non- 
local electrostatics based entirely on partial differential 
equations. 

In order to apply the equations to specific physical 
situations, the nonlocal dielectric function and thus the 
Green function of the differential operator C need to be 
determined. To keep the computations simple - since 
here we focus only on the basic conceptual features of our 
approach - we further assume that the dielectric function 
is isotropic, i.e., Q{r,r') « Q(r — r'), which is exact far 
from any boundary. The use of more general expressions 
is clearly permitted in our theory and is in fact needed for 
the treatment of realistic situations HI E 113. Within 
our approach, they lead to more complex Green func- 
tions and corresponding differential operators, and will 
be discussed in a detailed study later. 

A standard model for an isotropic nonlocal dielec- 
tric function is given by the so-called Fourier-Lorentzian 
model with a Yukawa-type kernel in real space. The cor- 
responding Green function reads 



1 e- 



r' = 



4tt r — r' 



(15) 



with the differential operator C being given by £ = A 
1/A^. With this choice ea. (fn|l reads 



(16) 



Note that due to Atps — 0, no differential operator ap- 
pears on the rhs of ea. l(TB|l . 

This result is interesting for two reasons. First, it illus- 
trates that in the limit A ^ 0, i.e. on length scales large 
compared to the scale of the orientational correlations, 
the local limit is recovered. Second, for the differential 
operator chosen, the form of the equation is apparently 
that of a Debye-Hiickel equation in which the role of the 
Debye-Hiickel screening length is played by the combi- 
nation \{ei/e^) and the potential ip plays the role of 
the density of mobile charges. We can thus interpret ip as 
a density of polarization charges in the bulk whose gradi- 
ent gives rise to the rotation free part of the displacement 
field D. 

We now turn to illustrate how our formulation of non- 
local electrostatics can be put to use. First, we consider 
the simplest case of a charge q placed at the center of a 



spherical shell of radius a. Inside the shell, we assume 
en = 1. This system serves as a model for ion solva- 
tion 9] . The equations for ip and can now be solved as 
follows. In the nonlocal case, the role of the Poisson equa- 
tion for (j) is taken over by the equation for tp according 
to eas. (|8l9|l . The -^-potential inside the shell is given by 

"00 = 4^7, with the same form given outside, 0s — 
The tangential boundary condition is trivially fulfilled, 
while the normal boundary condition at r = a leads to 
q = q' . Due to the radial symmetry of the problem, eq. 
(Unj reads 



£S 



T- (17) 



which is solved by 



^ ^(H-Aexp(-7r)) (18) 



47r£s£o r 



where the coefficients A and 7 follow from coefficient 
matching and the continuity of at the boundary. 



A: 



0n(a) 

£s — £1 sinh v 
£l V 



£s a 



£l A 



es 1 
£/ A 



(19) 



The electrostatic potential can be used to estimate the 
solvation energy of monovalent and divalent ions. From 
0Q,(r) and tpaif), with a € {fl, S}, we can easily compute 
the free energy of solvation for this setting as the differ- 
ence of the electrostatic energies in water and vacuum 
(where the local computation of course remains valid) 
from AG = ^ / {p(f> — p(f>vac}dr, where the integrals are 
split into integrals over and E. 

The result of this calculation is shown in Figure 2, 
where we have compared our results to a correspond- 
ing local computation (the Born- model). The correla- 
tion length A serves as an adjustment parameter; it is 
the only one in the theory. Like the ion radii, it can 
also be obtained from microscopic simulations. Our re- 
sult compares favourably to the experimental data taken 
from H3 ■ The value for A was taken to be 24.13A; for the 
ion radii we chose the values according to Aqvist [iJl ■ We 
also tested the set of Shannon radii 0, |23| without sig- 
nificant differences on our results. A detailed discussion 
of the choice of ion radii within nonlocal electrostatics 
can be found in |^ . 

The simple radially symmetric problem of ion sol- 
vation is, of course, not representative for the general 
character of the solutions of eg. ljlGI) . as the tangential 
boundary condition is trivially fulfilled in this case. In 
order to elucidate the effect of this boundary condi- 
tion, we consider the potential generated by a charge 
q placed at a distance d from a planar dielectric phase 
boundary, as sketched in Figure 3. Within local elec- 
trostatics, this problem can be solved by the method of 
image charges. The symmetry of the situation allows 
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FIG. 2: Solvation free energy of ions. 




FIG. 3: Geometry of the problem of a charge q placed at a 
distance d from a dielectric boundary. For a discussion, see 
the text. 



the use of cylindrical coordinates. At a point P inside 
fi, the local dielectric medium, the electrostatic poten- 
tial 4) is given by (t>n{p,z) = -j^f^^ + where 

Ri = ^ + (d - zf and i?2 = \/ + + are the 
lengths of the vectors pointing from the charge q located 
at z = d, and the image charge g' at z = — d to the 
point P. Inside the charge-free halfspace S the potential 
is given by (j>^ = T^^i^ ' ^ith the image charge q". 
The application of the boundary conditions for the nor- 
mal and tangential components then yields the relations 
between the image charges q+q' = ^q" , q—q' = q" so 
that the potential is determined explicitly in both halfs- 
paces. Inside ft, e.g., one has 



(t>nig,z) 



1 



ineoen \Ri £n + £s R2 



(20) 



In the nonlocal case, the eqs. for if) have the same form 
as those for (j) in the local case, so we assume the solu- 
tions to be similar. Under this assumption, the normal 
boundary condition remains unchanged while the tan- 



gential boundary condition leads to 

1 {q + q') 



do 



(21) 



at z = 0. This equation can be readily integrated along 
g and determines the potential at the dielectric bound- 
ary. For g, z —> 00, on the other hand, (/)x; must vanish. 
These two conditions then allow to compute the poten- 
tial 4>y: from ea. HKil) within E. Since this requires a nu- 
merical computation, we leave it for a future publication 
where we discuss the numerical treatment of our equa- 
tions. Here, we only give the lowest order effect nonlo- 
cality within E has on il. For g <C V^d, the electrostatic 
potential inside Q has the same form as in the local the- 
ory, but with a 'renormalized' bulk value of the dielectric 
function £s = esAoc — 2e(,{\^ /d^). Transverse variations 
of the permittivity in the nonlocal medium thus induce 
a change of the local permittivity proportional to (A/d)^ 
in the vicinity of the boundary, an effect which vanishes 
for A ^ 0, and for deeply buried charges, d — > 00. 

As a final remark on the applications of the nonlocal 
theory of electrostatics we mention the case in which mo- 
bile charges are present in the medium surrounding the 
cavity. Within the linear mean-field theory of local elec- 
trostatics, they can be described by Poisson-Boltzmann 
theory. This stays true within the nonlocal theory pre- 
sented here. The Boltzmann distribution of the charges 
simply modifies the rhs of eq. © . The nonlocal approach 
can then e.g. be used to quantify recent experimental re- 
sults of AFM measurements on force-deflection curves at 
charged mica substrates in water and solutions of mono- 
valent ions (23]. Within a simplified treatment of the 
dielectric function of water, the orientational effects gov- 
erned by A are found to be on the order of 10 nm [2^. We 
expect that a more realistic structural model for water, 
which will become computationally tractable due to our 
approach, will lead to a much improved description of 
water orientation near charged surfaces, at least in cases 
where the assumption of mean- field behaviour is justified. 

To conclude, we have presented a novel formulation of 
nonlocal electrostatics, which includes the effects of spa- 
tial dispersion in the dielectric permittivity on surfaces 
embedded in a solvent, by reformulating it in terms of a 
two-potential model. While the resulting equations still 
need to be solved numerically even for simple geometries, 
this task can now be performed by standard methods 
developed for partial differential equations. Due to the 
generality of eq. H14II , dielectric functions of greater com- 
plexity than the simple radially symmetric choice used 
here for illustrative purposes can be treated, provided 
they can be related to known Green functions. Work in 
this direction is under way. 
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